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General molecular dynamic approach, making possible direct calculation of eigen values and eigen 
functions for a quantum-mechanical system of an arbitrary symmetry is proposed. The method 
is based on analogy between discrete representation of the Schrodinger equation and the system 
of Newton equations describing dynamics of specially constructed virtual lattice. Few examples 
demonstrating the method capabilities are considered. 

PACS Numbers 73. 20. At, 73.20.Hb 

Molecular dynamics (MD) is a technique which is widely used for simulation vibration of atoms and molecules, 
growth of a crystal or interface, the interaction of an adatom and surface, and a wide range of other time-dependent 
phenomena. In the MD technique, the many-body classical equations of motion are solved as a function of time, and 
the physical process can be studied in real time if the instantaneous forces acting on an atom are given. The latter 
are obtained either using some devised potentials to mimic the complicated character of the forces in realistic objects 
jjj or from electronic-structure calculations (see, e.g. || and references therein). In the ab initio MD calculations the 
numerical integration of the Newton equations of motion for a given number of ions involves the forces derived, at 
each time step, from the instantaneous electronic configuration. 

While considerable efforts have been made to develop the quantum-mechanical electronic structure calculations 
and to use the obtained Hcllmann-Feynmann forces in the MD simulations (see, e.g. ||), little attention has yet been 
given to applying the MD technique in the quantum mechanics directly. Recently a fictitious Newtonian dynamics 
was introduced for electronic variables to make the ab-initio MD calculations more effective 0. The physics behind 
this method was analyzed in 0] . 

We present here a general MD approach, which allows direct calculation of eigen values and eigen functions of 
any stationary quantum-mechanical problem. To determine by ordinary methods the eigen functions are assumed 
to be expanded onto a complete set of orthogonal functions. These can be either atomic orbitals or plane waves 
depending on whether the localized or delocalized states arc described. Our method can be effectively used regardless 
the localization properties of the studied states, geometry and symmetry of a system. We restrict our consideration 
here by a one-particle multicenter problem although the proposed method can in principle be applied to many-particle 
problems as well. 

The method is based on discrete representation of the Schrodinger equation for the particle of mass m e in the 
potential U (x, y, z) 



H ( d 2 9(x,y,z) d 2 9(x,y,z) d 2 9(x.y,z) 



2m e \ dx 2 dy 2 
by a system of oscillator equations. Indeed 



+ (U(x,y,z)-E)-9(x,y,z) = 0, (1) 
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and analogously for the d ^g^'^ and d ^^qJi'^ 1 ■ Let the ^-function be defined on discrete space, i.e. x = n ■ 5 X , 
y = m ■ 5 y , z = I ■ 5 Z . Then 



*n,m,; = 9(n ■ 5 x ,m - 5 y ,l- 5 Z ), 
U n ,m,i = U(n ■8 x ,m- S y ,l ■ S z ) 



and Eq. ([!]) can be approximated as 
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' (2^n,m,/ ^n,m,/ + l ^'A,m.!-l) "1" U n ra i ' ^n,m.l 0; 

where Ki = ft 2 / (2m e 5 2 ) (i — x,y, z). Suggesting 

*n,m,j(*) = *n,m,i ' exp^v 7 !? • t) (3) 
we can rewrite (^|) and obtain the system of linear equations 

d2 ^'' (t) + K x ■ {29 n>m ,i(t) - 9 n+ i, m>l (t) - *«_i, m ,i(*)) + K y • (2*„, ro ,i(t) - *„, m+ i,i(0 ~ *„,m-i,i(*)) + (4) 

• (2*„, m ,i(t) - f n ,m,l+l(i) ~ $n,ro,l-lW) + C/^mJ • *n,m,i(<) = 

describing the harmonic vibrations of the unity-mass virtual particles in a lattice with nearest neighbor force constant 
Ki = h 2 /(2m e 5 2 ) and incite force constant U n . m ,i- The ^ n ,m.i{t) in this context means the displacement of the virtual 
particle in the cite (n,m,l) from its equilibrium position. Thus, from continues equation (|l|) for the one-particle W 
function we have arrived to the system of Newton equations (Q) describing the dynamics of the virtual lattice and the 
eigen functions of the quantum-mechanical problem are associated now with the eigen vectors of the virtual lattice. 
Therefore the described method one may call as the virtual lattice dynamics (VLD) method. 

For practical use of the VLD method one should construct first the virtual lattice and then follow the usual way of 
MD simulation of lattice dynamics. An important point concerns an a priori estimation of the accuracy of the VLD 
method since it results in the proper choice of spacings 5, . The accuracy can be clearly controlled using the definition 

Ej = (^S>j H V For the ground state Eq one may neglect the kinetic energy operator in the Hamiltonian keeping 

the potential term U(x, y, z) only. As the \I/ function is a smoother function than the potential one and its discreteness 
can be neglected, the accuracy can be estimated on a bases of approximation of the potential curve U(x, y, z) by the 
discrete function U n ,m,l- For 1-D case one may accept ^>j = Const over the ^>j function localization length Lj = Nj-S 
and then roughly estimate the relative accuracy of the VLD method as ~ where Nj is the number of the 

virtual particles inside the Lj . Similar arguments can also be applied to a 3-D case where relative accuracy roughly 
corresponds to the ratio ~ Nj /Nj , Nj and Nj being the numbers of virtual particles on the surface of the \&j 
function localization region and inside it respectively. Thus, the accuracy of the VLD method obviously increases 
with decrease of spacing 5 getting asymptotically exact. 

For illustration we consider the simple analytically solvable example from quantum mechanics textbook: single 
unity-mass particle in semi-infinite 1-D rectangular potential well (see Fig.l). The virtual lattice is represented by a 
simple linear chain of particles separated by 6 with interparticle coupling constant K = l/(2<5 2 ) (Ti = 1) and incite 
force constant U n = for n < 24 and U n = 0.1 otherwise. We have chosen the number of virtual particles inside the 
potential well N w = 23 (5 = 1) and hence the a priori estimated accuracy in the determination of eigen values must 
be about 4 per cent. At the t = the particles in the chain were randomly displaced and their subsequent vibrations 
were analyzed via Fourier transformation. According to (|J) the eigen energy Ej is given by wj, where u>j is the peak 
position in the Fourier spectrum ^ |^ n (o;)|. The Fourier spectrum ^2 \^n(^ 2 )\ shown in Fig.l reveals five fairly 

n n 

separated sharp peaks corresponding to the bound states within the potential well, and the continuum corresponding 
to extended states. The calculated eigen values for the bound states deviate from those obtained analytically by less 
than 2 per cent giving evidence for quite reasonable a priori estimation of the accuracy. In case of N w = 48 (<5 = 1/2) 
the deviation decreases down to ~1 per cent and obviously will tend to zero if 5 — > 0. 

After determination of the eigen values one can calculate the corresponding eigen functions. In nondegenerate case 
the eigen function $ j (x) can be approximated as 

y(x)~V{n-6)=C- 1 -Re[y n (uj j )}, (5) 

N 2 

where C = ^2 [^ni^j)]} is normalization factor, TV is the total number of virtual particles in the chain. Another 

n=l 

and more exact though more time consuming way to determine the ^j(x) is as follows: a virtual particle at a 
nonsymmetric position m is initially excited by a harmonic force with eigen frequency uij = \fE~j. To reach the 

steady state condition of the excitation the small phenomenological damping 7 • was introduced into the 

motion equations (Q). After few tens of vibrations the spatial character of the virtual particles vibration around 
in — th particle is recorded and the excitation force field is made identical to the particles displacement field. The 
procedure is repeated until no noticeable difference is observed between the excitation and the particles displacement 
fields. The final spatial character of the excitation must be identical to the corresponding \&j(ar) function since no 
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other states can be excited if there is no degeneracy. Thus obtained ^j(x) functions (see, e.g. Fig.l) are in perfect 
agreement with the theoretical ones. 

To check if the eigen state is degenerated a different virtual particle have to be chosen for the initial excitation 
with the same frequency u)j. The state is nondegenerated if the new determined &j(x) function coincides with that 
determined previously for the m — th virtual particle as the initially excited one. Otherwise one should determine the 
degree of degeneracy by taking different virtual particles for the initial excitation and then construct the orthogonal 
combinations from all determined ^ j (x) functions for the given uij . 

To illustrate the applicability of the VLD method to a more complicated system we calculated the energy levels and 
the charge density distribution in the bulk and near the surface of the 2-D two-atomic crystal band with impurities 
(Fig. 2). As an ordinary MD technique the VLD method is restricted by the number of particles under treatment. 
Being able to handle effectively about a billion of virtual particles we describe the crystal band consisting of 10 x 40 
potential wells (atoms) with an a priori estimated accuracy of few per cent. Virtual lattice in this case is a simple 
quadratic one extended outside the crystal band in the x-direction. Cyclic boundary conditions in both x- and y- 
directions for the virtual lattice vibrations were applied. Panels (a) and (b) in Fig. 3 show the shape of the atomic 
potential. The potential parameters were chosen in a way to form a surface electronic band near the middle of the 
band gap in the energy spectrum. The latter obtained via Fourier analysis of the virtual lattice vibrations is shown 
in Fig. 3c. The potential of impurity atoms (denoted by stars in Fig. 3) was chosen to be only slightly different from 
that of one of the host atoms. 

Using the excitation in the form of standing wave with different period the dispersion of the electronic excitations 
over the Brillouin zone was obtained (Fig. 4). Different dispersion curves in Fig. 4 can be associated with electronic 
bands of different symmetry, what is not a subject of our paper, however. We are interested basically in discrimination 
of the states with regard to the wave vector k y . 

Spatial charge distribution Pj(x,y) ~ ^j(x,y) in different eigen states corresponding to the Brillouin zone centre 
k y = are shown in Fig. 5 for both pure and impurity crystals. The potential defect related to impurity atom could 
not result in the localization in the bulk of the crystal, but turned out to be big enough to form the localized states 
at the surface (see Fig. 6). The latter are strongly enhanced around the impurity clusters near the surface. The VLD 
method allows also to study an influence of surface roughness on electronic states. Charge distribution in the bulk 
eigen states disturbed by the surface roughness is shown in Fig. 7. The roughness induced localized state at the surface 
is presented in Fig. 8. 

The considered example showed that the VLD method can be effectively used for numerical studies of one-particle 
eigen states in inhomogeneous and low symmetry systems. The method seems quite adequate for direct investigation 
of defect and/or impurity states in compensated semiconductors when the interaction between one-particle localized 
states of donors can be neglected. In case of the nonvanishing interaction the VLD method can still be applied via 
proper iteration procedure. Worth to note one important feature of the VLD method: remaining asymptotically exact 
it can be generalized on the case of few particles in multicenter potential allowing to solve excitonic problem in a 
cluster of an arbitrary symmetry. 
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Figure captures 

Fig.l. One-particle states in the semi- infinite rectangular potential well (shown by thick solid lines in the right panel) 
U n = if n < 38 and U n = 0.2 otherwise. The energy levels correspond to the peak positions in the Fourier 
spectrum ^ |^ , „(w 2 )| (right panel). The eigen functions calculated for three lowest lying eigen states are shown 
in the left panel. 

Fig. 2. Two atomic 2-D crystal with impurities. Axes labels are given in <5 units thus corresponding to the virtual 
particle number. Large and small black points denote different potential wells corresponding to two different 
types of the host atoms. The atoms with deeper potential are randomly substituted by impurity atoms (white 
points). Atomic potential profile is shown in Fig. 3. 

Fig. 3. Potential profile of the crystal shown in Fig. 2: a) along x-axes at y=6; b) along y-axes at x=133; c) Fourier 
spectrum of the virtual lattice vibration showing the valence (V) and conduction (C) bands in the bulk, surface 
band (S), and localized states (L) due to impurity atoms at the surface. 

Fig. 4. Dispersion of eigen energies over the Brillouin zone for the 2-D two atomic crystal a) pure; b) with impurities 
(Fig. 2); c) pure crystal with rough surface (Fig. 7a). Surface band is labelled by S. Open circles denote eigen 
energies for which the cigen functions were calculated (Figs 5-8). 

Fig. 5. Charge density distribution pj(x, y) <~ ^(x, y) calculated for some eigen states of the pure ((a), (b), and (c)) 
and the impurity ((d), (e), and (f)) crystals: (a) and (d) correspond to the state 1, (b) and (e) - 2, (c) and 
(f) -3, denoted in Fig. 4. The states (a), (b), (d), and (e) are obviously the bulk ones while (c) and (f) are the 
surface ones. 

Fig. 6. Charge density distribution pj(x,y) <~ ^ 2 -{x,y) in the localized surface state (right panel) of the impurity 
crystal (labelled as IL in Fig. 4) in comparisom with the extended surface state (left panel) of the pure crystal 
(3 in Fig. 4). Open circles in the right panel denote impurity atoms. 

Fig. 7. ^(x,y) calculated for the pure crystal with rough surface (a). Panels (b), (c), (d) correspond respectively to 
the states 1, 2, 3 in Fig. 4. 

Fig. 8. ^"j(x,y) calculated for the roughness induced localized (b) surface state (RL in Fig. 4) in comparison with the 
extended surface state (a) of the pure crystal (4 in Fig. 4). 
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